
#include "../common/improve_div.h"
#include "../common/matrix_op.h"
#include "../common/mem.h"
#include "../common/error.h"

#include <stdlib.h>
#include <memory.h>
#include <float.h>

int max_unmoved_element(const int* unmove_vector,
						const int unmoved_count,
						const elem* vector,
						const int vector_size)
{
	int i;
	elem max_value;
	int max_index;
	int unmoved_index;

	if (vector_size <= 0)
		return -1;

	max_index = -1;
	max_value = -DBL_MAX;

	for (unmoved_index = 0; unmoved_index < unmoved_count; unmoved_index++)
	{
		i = unmove_vector[unmoved_index];
		/*  Unmoved and larger than the maximum value? */
		if (vector[i] >= max_value)
		{
			max_index = unmoved_index;
			max_value = vector[i];
		}
	}

	return max_index;
}

void calc_scores(const ModularityMatrix* Bg, const int group_size,
						 const int* s, elem* scores)
{
	int i;
	Vector* x;
	Vector* temp;
	elem score;
	const elem* Kg = Bg->Kg;
	const elem M = Bg->M;

	x = allocate_Vector(group_size);
	temp = allocate_Vector(group_size);

	/*  temp = S */
	for (i = 0; i < group_size; i++)
		temp->data[i] = (elem) s[i];

	/*  x = Bg * temp */
	modularity_mat_mlpl_by_vector(Bg, temp, x, 0);
	/* matrix_vector_mlpl(Bg, temp, x); */

	/*  Calculate scores, and their maximum */
	for (i = 0; i < group_size; i++)
	{
		/*  TODO: Reverse S[i] sign? */
		score = (-2) * (s[i] * x->data[i] + Kg[i] * Kg[i] / M);
		scores[i] = score;
	}

	free_vector(x);
	free_vector(temp);
}

void update_unmoved_scores(const ModularityMatrix* Bg,
						   const int* unmoved_vector, const int unmoved_count,
						   int* s, int moved_index /* K */, elem* scores)
{
	int i;
	int k;
	int unmoved_index;
	elem partial_computation;

	k = moved_index;
	scores[k] *= -1.0;
	partial_computation = 4 * s[k];

	/*  Note:  */
	/*  In the loops below (k,i) is the same as (i,k) */
	/*  because of symmetry, but more efficient to get. */

	for (unmoved_index = 0; unmoved_index < unmoved_count; unmoved_index++)
	{
		i = unmoved_vector[unmoved_index];
		scores[i] -= partial_computation * s[i] * 
			modularity_mat_get_cell(Bg, i, k, 0);

		/* if (modularity_mat_get_cell(Bg, i, k, 0) != modularity_mat_get_cell(Bg, k, i, 0))
						fatal_error("BAD!");*/
	}

}

elem improve_div_iteration(const ModularityMatrix* Bg, const int group_size,
						   int* s, elem* scores, int* unmoved_vector,
						   elem* improvement)
{
	elem result;
	int max_score_index;
	int i,j;
	int max_improve_index = -1;
	elem max_improve_value;
	int* indices;
	int unmoved_count;
	int max_unmoved_score_index;

	indices = (int*) mem_alloc_and_zero(group_size,sizeof(int));
	result = 0;
	max_improve_value = -DBL_MAX;

	/*  Populate the "unmoved" vector */
	for(i = 0; i < group_size; i++)
		unmoved_vector[i] = i;

	/*  We start by all indexes being "unmoved" */
	unmoved_count = group_size;

	/*  Compute Delta-Q for each move */
	calc_scores(Bg, group_size, s, scores);

	/*  FOR i = 0 to (Ng - 1) do ... */
	for (i = 0; i < group_size; i++)
	{
		/*  NOTE: Scores are "updated" instead of "calculated" at the end of the */
		/*  loop body, as specified by the project's description for efficient */
		/*  score computation. */
		/*  */
		/*  First scores were calculated before this function. */

		/*  "max_score_index" == "j'" in the documentation */
		max_unmoved_score_index = max_unmoved_element(unmoved_vector, unmoved_count,
			scores, group_size);

		max_score_index = unmoved_vector[max_unmoved_score_index];

		/*  Move vertex with maximal score */
		indices[i] = max_score_index;
		s[max_score_index] *= -1;

		if (i == 0)
			improvement[0] = scores[max_score_index];
		else
			improvement[i] = improvement[i - 1] + scores[max_score_index];

		/*  Update moved vector (remove "moved" index) */
		unmoved_count--;
		unmoved_vector[max_unmoved_score_index] = unmoved_vector[unmoved_count];

		/*  Update the scores instead of re-calculating them (more efficient) */
		update_unmoved_scores(Bg, unmoved_vector, unmoved_count, s,
			max_score_index, scores);
	}

	/*  Find maximal improvement in S and update S accordingly */
	
	/*  Get max improvement */
	max_improve_index = max_element(improvement, group_size);

	debug_output("Max improve index: %i\n", max_improve_index);

	for (i = (group_size - 1); i >= (max_improve_index + 1); i--)
	{
		j = indices[i];
		s[j] *= -1;
	}

	if (max_improve_index == (group_size - 1))
		result = 0;
	else
		result = improvement[max_improve_index];
	
	mem_free(indices);
	return result;
}

elem improve_div(const ModularityMatrix* Bg, int* s)
{
	elem delta_Q;
	elem total_delta_Q;
	elem* scores;
	int* unmoved_vector;
	elem* improvement;
	int group_size = Bg->Ag->n;

	unmoved_vector = mem_alloc(group_size, sizeof(int));
	scores = mem_alloc(group_size, sizeof(elem));
	improvement = mem_alloc(group_size, sizeof(elem));

	total_delta_Q = 0;

	do
	{
		delta_Q = improve_div_iteration(Bg, group_size, s, scores,
			unmoved_vector, improvement);

		if (IS_POSITIVE(delta_Q))
			total_delta_Q += delta_Q;

	} while IS_POSITIVE(delta_Q);

	mem_free(scores);
	mem_free(improvement);
	mem_free(unmoved_vector);

	return total_delta_Q;
}

